#include <float.h>
#include <math.h>

#define SQRT_3_4PI M_2_SQRTPI*9/32

typedef struct {
    float x;
    float y;
    float z;
    float i[3]; //voxel's light intensity - RGB
} Voxel;

typedef struct {
    Voxel origin;
    float i[3];
} LightSource;

typedef struct {
    Voxel vertex[8];
} BoundingBox;

Voxel newVoxel(float x, float y, float z) {
    Voxel v;
    v.x = x;
    v.y = y;
    v.z = z;
    return v;
}

int isInside(Voxel v, BoundingBox b) {
    int i;
    float min[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
    float max[3] = {FLT_MIN, FLT_MIN, FLT_MIN};
    //Find minimum and maximum bound box coordinates
    for(i=0;i<8;i++) {
        //Minimum x
        if(b.vertex[i].x < min[0]) {
            min[0] = b.vertex[i].x;
        }
        //Maximum x
        if(b.vertex[i].x > max[0]) {
            max[0] = b.vertex[i].x;
        }
        //Minimum y
        if(b.vertex[i].y < min[1]) {
            min[1] = b.vertex[i].x;
        }
        //Maximum y
        if(b.vertex[i].y > max[1]) {
            max[1] = b.vertex[i].x;
        }
        //Minimum z
        if(b.vertex[i].z < min[2]) {
            min[2] = b.vertex[i].x;
        }
        //Maximum z
        if(b.vertex[i].z > max[2]) {
            max[2] = b.vertex[i].x;
        }
    }
    //Test if given voxel's coordinates are included in coordinate ranges found
    int insideX = min[0] <= v.x && v.x <= max[0];
    int insideY = min[1] <= v.y && v.y <= max[1];
    int insideZ = min[2] <= v.z && v.z <= max[2];
    return (insideX && insideY && insideZ);
}

float distance(Voxel v1, Voxel v2) {
    return sqrt(pow(v1.x - v2.x, 2) + pow(v1.y - v2.y, 2) + pow(v1.z - v2.z, 2));
}

//From this point onward, I have no idea what I'm doing
float p(float coord, float radius) {
    return SQRT_3_4PI * coord/radius;
}

float SH(float r, Voxel v) {
    float s = 0.25 * M_2_SQRTPI;
    float p_x = p(v.x, r);
    float p_y = p(v.y, r);
    float p_z = p(v.z, r);
    
}

float s(int j) {
    return j;
}

float f(float s_j) {
    return s_j;
}

float yml(float s_j) {
    return 2*s_j;
}

float MC(int N) {
    int j;
    float s_j, S=0;
    for (j=1;j<=N;j++) {
        s_j = s(j);
        S += f(s_j)*yml(s_j);
    }
    return 4*M_PI*S/N;
}